Using Correlation Integrals to Characterize 3D Stellar Orbits 
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ABSTRACT 

In an effort to more fully understand the variety of stellar distribution functions that can be 
used to construct models of realistic galaxies, the correlation integral method for orbit character- 
ization that previously has been introduced by Grassberger & Procaccia (1983a) and Carnevali & 
Santangelo (1984) is examined in considerable detail. The broad utility of the method is validated 
and demonstrated by applying it to orbits in a number of different, previously studied test cases 
(ID, 2D, and 3D; nonrotating and rotating). At the same time, the correlation integral method 
is compared and contrasted with other more traditional characterization tools such as Lyapunov 
exponents and surfaces of section. The method is then extended to orbits in a previously un- 
examined rotating, 3D bar potential. The correlation integral method is found to be a simple 
and reliable way to quantitatively categorize orbits in virtually any potential. It is recommended 
that it be broadly adopted as a tool for characterizing the properties of orbits and, by extension, 
stellar distribution functions, in all Hamiltonian dynamical systems. 
Introduction 



1.1. Background 

The study of orbits in three-dimensional (3D) 
gravitational potentials is of broad astrophysical 
interest. Whether the orbits under investigation 
are those of asteroids in the solar system (e.g., 
Pilat-Lohinger et al. 1999), stars in a globular 
cluster or a galaxy (e.g., Carpintero et al. 1999), or 
galaxies in a cosmological simulation (e.g., Colpi 
et al. 1999), methods of characterizing the orbits 
can be extremely useful, especially when an at- 
tempt is made to understand what the "typical" 
behavior is of a very large collection of orbits in a 
particular physical system. These methods, some 
of which are discussed below in §3, are most useful 
when they are quantitative, reliable, and applica- 
ble to a wide variety of systems. 

Here, the analysis will focus on stellar orbits in 
models of time-invariant galactic potentials. With 
this in mind, a brief review of some basics of stel- 
lar dynamics is presented. Stars in galaxies are 
presumed to form a coUisionless system, that is, 
the forces felt by any one star are due only to the 
mean forces produced by all other stars. The dis- 
tribution function, /, of stars must therefore sat- 



isfy the coUisionless Boltzmann equation (Binney 
& Tremaine 1987), 
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where the Wi are the phase space coordinates, e.g., 
{x,y,z,x,y,z). A by-product of this equation is 
that any time-independent function of phase space 
coordinates, known as an integral of motion, must 
be a valid solution of eq. (1). This fact is incor- 
porated into the Jeans Theorem: any steady-state 
solution of eq. (1) must be a function of integrals 
of motion and any function of integrals of motion 
must likewise solve eq. (1) (Binney & Tremaine 
1987). A stronger, qualified statement can be 
made. If all orbits in a given potential are reg- 
ular^ then any solution to eq. (1) is a function 
of three integrals of motion. This is known as 



^We use the term "regular" to describe orbits that respect 
a number of isolating integrals of motion that is greater 
than or equal to the degrees of freedom of the orbit (Bin- 
ney & Tremaine 1987). For example, an orbit in a 2D 
potential must have at least two isolating integrals in or- 
der to be considered regular. On the other hand, we use 
"quasi-ergodic" as a blanket term for both stochastic and 
semistochastic orbits; irregular orbits are fully ergodic. Er- 
godic and quasi-ergodic orbits respect only one complete 



1 



the Strong Jeans Theorem (Binney & Tremainc 
1987). Indeed, in potentials that have analytically 
prescriptible integrals of motion, such as spherical 
or axisymmctric systems, a variety of distribution 
functions have been investigated (e.g.. King 1966; 
Kalnajs 1976). 

However, most realistic models of galaxies are 
not as simple as those discussed above, and differ- 
ent tactics must be used to investigate the prop- 
erties of the stellar distribution functions that are 
associated with these more realistic systems. And 
while an in-depth discussion of the creation of 
such models is beyond the scope of this paper, a 
short overview is given here. N-body simulations 
(e.g.. Miller & Smith 1979; Hohl & Zang 1979; 
Combes & Sanders 1981; Miller et al. 1982; Pfen- 
niger & Friedli 1991) have been widely used to 
generate stellar distribution functions for steady- 
state galaxy models. By design, many of these 
simulations have begun with axisymmetric dis- 
tributions of point particles, then the dynamical 
systems quickly deform to bar-shaped objects. 
Most of these studies have focused on the global 
evolution of the bar, but some (Miller & Smith 
1979; Pfenniger & Friedli 1991) have also looked 
at constituent orbits of the bars and have found 
evidence for nonanalytical integrals of motion. 
Performing numerical integrations of individual 
orbits in analytical potentials provides further ev- 
idence for these nonanalytical integrals in studies 
designed to find periodic orbits (e.g., Hcislcr et 
al. 1982; Magnenat 1982; Mulder & Hooimeyer 
1984; Pfenniger 1984; Martinet & de Zeeuw 1988; 
Hasan et al. 1993). Uncovering periodic orbits 
has been very important in the context of stud- 
ies of 3D galaxy models because Schwarzschild 
(1979, 1982) developed procedures for creating 
self-consistent potential-density pairs (and associ- 
ated steady-state distribution functions) once the 
periodic orbits for a given 3D potential arc known. 
These techniques have been extended to model re- 
alistic stellar dynamical systems (e.g., Rix et al. 
1997: van dcr Marcl et al. 1998; Cretton et al. 
1999; Cretton, Rix, & de Zeeuw 2000). 

While the importance of regular orbits has 
been the focus of many stellar dynamics re- 
search projects, the importance and likely exis- 
tence of quasi-ergodic orbits in realistic galaxy 

isolating integral of motion, either the specific energy, e, or 
the Jacobi constant, ej. 



potentials also has been discussed (Goodman & 
Schwarzschild 1981; Binney 1982b; Habib et al. 
1997; Vahuri & Merritt 1998). While the previ- 
ously mentioned studies involved analytical poten- 
tials, Barnes & Tohline (2001) have also found that 
quasi-ergodic stellar orbits are quite common in 
a numerically-created, gaseous bar. Quasi-ergodic 
orbits also appear in N-body simulations of galac- 
tic bars (see, for example, Sparke & Sellwood 
1987). A specific example of a situation in which 
quasi-ergodic orbits play an important role is in 
galaxy models with massive central objects (e.g., 
Udry & Pfenniger 1988, Valluri & Merritt 1998). 
In a detailed analysis of the influence that mas- 
sive central objects have on stellar orbits, Gerhard 
& Binney (1985) suggest that sufficiently massive 
central objects will scatter into quasi-ergodic or- 
bits trajectories that initially: 1) are regular; 2) 
pass near the central object; and 3) support non- 
axisymmetric distributions, such as bars. With 
the current interest in understanding active galac- 
tic nuclei as well as the recent discovery of a link 
between the masses of central objects and central 
stellar velocity dispersions (Ferrarese & Merritt 
2000; Gebhardt et al. 2000a), more accurate mod- 
eling of the distribution functions of such systems 
(along the lines of van der Marel et al. 1998) 
could provide some insight into the observational 
evidence. 

In an effort to better characterize orbits (both 

regular and quasi-ergodic) in astrophysically in- 
teresting potentials, the relative utility of the cor- 
relation integral method has been examined. As 
pointed out by Grassbcrgcr & Procaccia (1983a) 
and Carnevali & Santangelo (1984), the correla- 
tion integral method is a flexible and accurate 
characterization technique that has been widely 
utilized by physicists but has been virtually ig- 
nored by the stellar dynamical community. This 
technique distinguishes orbits based on the num- 
ber of isolating integrals of motion that are re- 
spected by any given orbit. This makes it a useful, 
quantitative tool for examining stellar distribution 
functions that would arise in potentials such as the 
ones mentioned above. Regular and nonregular 
orbits can be differentiated. Within the regular 
class, periodic and quasi-periodic orbits can also 
be distinguished from one another. Additionally, 
this technique can distinguish between 3D orbits 
that respect two isolating integrals and those that 
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respect only one. This ability to distinguish in 
a clear and quantitative fashion between vari- 
ous types of orbits makes the correlation integral 
method particularly well suited to involvement in 
studies of 3D potentials where quasi-ergodic or- 
bits are likely to play an important role. The 
primary objective of this paper is to demonstrate 
the quantitative reliability and broad utility of 
the correlation integral method by applying it as 
a characterization tool to orbits in a wide variety 
of potentials - most of which are familiar to the 
stellar dynamics community - and by comparing 
and contrasting it to other techniques that have 
been broadly used to characterize orbits. This will 
place the community in a position to effectively 
utilize this tool in a wide variety of problems. 

The remainder of this paper is divided as fol- 
lows. The different models (one map and several 
potentials) arc introduced individually in §2. Sec- 
tion 3 contains summaries of the numerical inte- 
gration technique and various orbital characteriza- 
tion methods that arc utilized in this study. Tests 
of and results acquired with the correlation inte- 
gral method are presented in §4. Also discussed 
in this section arc the results from a previously 
unexamined 3D potential. Finally, §5 presents a 
summary of the findings. 

2. Models 

In order to illustrate the basic utility of the cor- 
relation integral method, as well as the accuracy 
of this particular implementation of the method, 
several systems that previously have been well 
studied and characterized using other techniques 
are presented. 

2.1. Henon Map 

The simplest system to attack using the correla- 
tion integral method is one that can be contained 
in a 2D phase space. The Henon map is one such 
system that exhibits nontrivial behavior. It is pre- 
scribed by the following set of iterative equations 
(Henon 1976): 

Xi+i =yi + l- axi, (2) 

and 

yi+i = bxi, (3) 



where a and b arc constants. Rather than solv- 
ing a pair of coupled differential equations for 
position and velocity, phase space is occupied by 
choosing initial x and y values, then iterating back 
and forth between these two algebraic expressions. 
When a = 1.4 and b = 0.3, the Henon map has 
a strange attractor in phase space, that is, the 
dimensionality of occupied phase space is nonin- 
teger (fractal). 

2.2. Richstone Potential 

Following the lead of Carnevali & Santangelo 

(1984), orbits in the scale- free, logarithmic po- 
tential referred to here as the Richstone potential 
(Richstone 1980, 1982) are also studied. (A tri- 
axial version of this potential has been studied in 
Binney 1981 and is known as Binney's potential.) 
One nice feature of this potential is that it can be 
used to study either 2D or 3D orbits (4D and 6D 
phase spaces, respectively). To investigate fully 
3D orbits, the following form of the potential is 
utilized, 

$r(x, y,z) = f In (^x' + y^ + ^+ , (4) 

where vq is the constant circular speed for the po- 
tential, g is a measure of the flattening of the po- 
tential, and Rc is a core radius. There are three 
possible potential shapes: 

• q> 1 - the potential is prolate spheroidal; 

• q = l - the potential is spherical; 

• q <1 - the potential is oblate spheroidal. 

As in Richstone (1982), the parameters are set as 
q = 0.75 and Rc = 0.1. 

In order to examine 2D orbital motion in the 
Richstone potential, thereby reducing the analysis 
from a 6D to a 4D phase space, first transform 
eq.(4) to cyhndrical coordinates {R,(j),z). It is 
clear that $r is axisymmctric, so the Lagrangian 
is cyclic in </>, which means that the ^;-component 
of angular momentum is constant. With this 
restriction, orbits in the 2D Richstone potential 
move under the influence of the following effective 
potential, 

$R,eff(i?,.) = |ln(i?^ + ^+i?^)+^. (5) 
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It is this 2D potential that Carncvah & Santan- 
gelo (1984) used to validate their implementation 
of the correlation integral method. 



2.3. Henon-Heiles Potential 

The Henon-Heiles potential (Henon & Heiles 
1964) is a well-known potential that supports both 
regular and quasi-ergodic orbits. This potential 
has the form, 

*HH(i?,2) = ^(i?' + 22 + 2i?2^-^23). (6) 

For energies e < 0.01, almost all orbits in this 
potential are regular. As the energy is increased, 
more and more of phase space is occupied by quasi- 
ergodic orbits. At an energy e = 1/6, the situation 
is reversed in that most orbits are quasi-ergodic. 



2.4. Stackel Potential 

Since Stackel potentials are separable, all orbits 
in such potentials are regular. This feature makes 
orbits in Stackel potentials ideal subjects for test- 
ing characterization techniques. The Stackel po- 
tential that is used in this study has the form (de 
Zeeuw & Lynden-Bell 1985), 



l + (f)^ + (f)^ + (f)^ 



(7) 



where vq = 1.0, a = 1.0, b — 0.8, and c — 0.6. 
This is a triaxial potential where the minor axis 
lies along the z-axis, the intermediate axis lies 
along the t/-axis, and the major axis lies along the 
a;- axis. 



2.5. Cazes Bar Potential 

In an effort to investigate the fission hypoth- 
esis for binary star formation, Cazes (1999) has 
produced two models of rapidly rotating, steady- 
state, triaxial, gaseous density distributions using 
a numerical hydrodynamics code. These models 
have compressible (n = 3/2) polytropic equations 
of state and nontrivial internal flows, including 
two relatively weak standing shocks. A detailed 
description of these models appears in Cazes & 
Tohline (2000). The hydrodynamics simulation 



that created these models was performed with di- 
mensionless units. This allows the results to be 
scaled to a variety of systems simply by choosing 
an appropriate mass and length scale (Williams 
&: Tohline 1987). For example, choosing a mass of 
1.0 Mq and a length scale i?eq ~ 4.57 AU describes 
a protostellar object of density p w 10~^gcm~^ 
(Cazes 1999). Alternatively, for a mass of 10^° Mq 
and length scale Req « 2kpc, the now galactic- 
sized object has a density p w 10~^^ g cm""^. With 
this galactic scaling in mind, the Cazes & Tohline 
(2000) "Model B" is adopted as an example of 
a nontrivial, rotating, bar-like structure with a 
self-consistent and realistic potential-density pair. 
Hereafter, this model will be referred to as the 
"Cazes bar." 

Since the Cazes bar was created numerically, 
values of various quantities such as mass density 
and the gravitational potential are specified at dis- 
crete locations on a computational grid; specifi- 
cally, Cazes & TohUne (2000) used a 128 x 64 x 256 
cylindrical (i?, z, (j)) grid. For reasons that will be 
made clear in a later section, it is difficult to accu- 
rately evaluate Lyapunov exponents and the cor- 
relation integral on such a coarse, discrete grid. In 
analyzing orbits in the Cazes bar potential, an an- 
alytical potential constructed along the lines dis- 
cussed in Barnes & Tohline (2001) to closely re- 
semble the numerical Cazes bar potential will be 
adopted instead. Specifically, the following ana- 
lytical potential will be used: 



^aCB(a;,2/, z) 



Nil- I 



L2 



+ 



f y 



\qR 



L2 



QzZlu 



(8) 



where iV is a normalization factor; q determines 
the strength of the bar-like distortion in the equa- 
torial plane; Qz determines the strength of the bar 
distortion in the x — z plane; and a, /3, and 7 
are exponents whose values are to be determined. 
Rl2 is the distance along the major axis from the 
center of the bar to its L2 Lagrange point, and 
Zlim IS a vertical scale height. Only values of q, 
Qz, Ziim, and i?L2 for which the a;-axis coincides 
with the major axis of the bar will be considered; 
the y-axis is then the intermediate axis; and the 
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z-axis is along the bar's minor axis. The angular 
velocity of the bar and the value of the potential 
at the center are taken from the numerical Gazes 
bar, that is n = 0.5218 and $min = -1.018, re- 
spectively. When z is set to zero, this potential 
reduces exactly to the 2D analytical potential that 
was studied in Barnes & Tohlinc (2001). 

In an effort to illustrate how well the analytical 
approximation to the Gazes bar potential $aCB 
matches the potential that was derived numer- 
ically in the Gazes & Tohline (2000) hydrody- 
namical simulation. Fig. 1 presents equipotential 
contours from both the numerical Gazes bar and 
the analytical Gazes bar. Figure la shows a slice 
of the numerical Gazes bar along the positive half 
of the major axis in the meridional plane. The cor- 
responding slice from the analytical Gazes bar is 
shown in Fig. lb. The parameters for the analyt- 
ical potential pictured here are: N = 0.7, q = 0.8, 
q, = 1.5, i?L2 = 1.36, zii^ = 0.65, a = 4, /3 = 4, 
and 7 = 1.7. Equipotential contours in the y — z 
plane of the numerical Gazes bar are illustrated in 
Fig. Ic. Figure Id displays the corresponding an- 
alytical Gazes bar plot. Figures le and If display 
equipotential contours in the equatorial planes of 
the numerical Gazes bar and analytical Gazes bar, 
respectively. The properties of this equatorial- 
plane, 2D potential were investigated thoroughly 
in Barnes & Tohline (2001); it will be used be- 
low as one of the test potentials. As a check on 
the applicability of the analytical Gazes bar as a 
substitute for the numerical Gazes bar, several in- 
tegrations with identical initial conditions in both 
potentials have been performed. In every case, the 
resulting orbital projections had similar (although 
not exact) morphologies. 

3. Integration and Characterization Tools 
3.1. Numerical Orbit Integration 

Over the past few decades, numerical integra- 
tions of the equations of motion, have become the 
standard way to investigate stellar orbits in galax- 
ies. In general, the equations of motion have the 
following vector form, 

f = -V$ - [fi X (O X x)] - 20 X f , (9) 

where O is the angular velocity of a rotating frame 
of reference, x is the position vector in that frame, 



and is the gravitational potential of the system. 
When a rotating potential is investigated in this 
paper, the angular velocity vector always points 
in the z-dircction. 

Orbits are calculated with a Verlet integration 
scheme (Verlet 1967). This is straightforward to 
implement for non-rotating potentials, such as the 
2D and 3D Richstone potentials. When a rotat- 
ing potential is examined, the Verlet scheme is 
modified so that two Verlet steps are performed 
per fixed timestep. This is done because Gorio- 
lis terms must be included in the accelerations. 
In order to provide optimum values of velocities 
for evaluation of the Goriolis terms, a first Verlet 
step is used to obtain a first estimate of the ve- 
locities, but particle positions and velocities are 
not permanently updated at this step. Then, for 
the same timestep, the second Verlet step returns 
and updates a more accurate subsequent position 
and velocity. This procedure provides an ade- 
quate level of energy conservation (Ae/e < 10~^ 
over lO'' timcstcps) for the purposes of this inves- 
tigation. The need for this level of energy con- 
servation is the reason that an analytical formula 
approximating the Gazes bar potential has been 
adopted in preference to the original, numerically 
prescribed $cb- Using an interpolation scheme 
with the numerical Gazes bar, the energy error 
(with the same timestep) is Ae/e 10"^. While 
this is adequate for a determination of the overall 
shapes of various orbits, the correlation integral 
method requires that the orbit integration scheme 
provide much better energy conservation. In prin- 
ciple, the numerical potential could be utilized if 
a smaller integration timestep were used. In prac- 
tice, though, this made total integration times too 
long to be practical. This retreat to an analyt- 
ical form of the Gazes bar potential should not 
be taken as a withdrawal from the assertion that 
the correlation integral method is applicable to 
numerically derived potentials. On a faster ma- 
chine (or with an improved interpolation scheme) 
the computational costs would not have posed as 
great a problem and a numerical potential could 
have been used. 

Initial conditions for the orbits presented in this 
paper have not necessarily been chosen to repre- 
sent a complete sample of phase space. Instead, 
the orbits discussed here illustrate different (but 
not all) families that exist in the various poten- 
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tials. For example, in the Hcnon-Heiles poten- 
tial, two regular and two quasi-ergodic orbits are 
followed to illustrate their respective differences. 
Orbits in the 3D analytical Cazcs bar potential 
are populated under the Restriction Hypothesis 
of Barnes & Tohline (2001). The idea behind the 
Restriction Hypothesis is that if stars form from 
gas in a system, the initial velocities of the stars 
will be determined by the gas flow at the points 
of formation. In practice, the Restriction Hypoth- 
esis is imposed by choosing initial positions that 
coincide with cylindrical grid points where the gas 
velocities (ux.gas, I'y.gas, ^^z.gas) from the numerical 
Gazes bar are known and most definitely does not 
represent a complete sampling of phase space. 

3.2. Orbit Morphology 

Perhaps the simplest way to categorize any 
orbit is according to its shape in configuration 
{x, y, z) space. Examples in 2D potentials are the 
"banana" orbits discussed in Binney (1982a), the 
4/1 orbits described in Contopoulos (1988), and 
the "bow tie" orbits presented in Barnes & Tohline 
(2001). For 3D orbits, projections of an orbit onto 
the three principal planes of the chosen configu- 
ration space coordinate system can illustrate the 
orbital morphology and therefore also can be use- 
ful when categorizing orbits. For example, if the 
X — y plane projection of an orbit appears to have 
definite circulation about the center (i.e., the z- 
component of the orbital angular momentum, Lz, 
is never zero), and the y — z and x — z projections 
appear as rectangular areas, the orbit is called a 
2-axis tube orbit (Binney & Trcmaine 1987). The 
benefit of characterizing orbits in this way is that 
such a descriptive name confers a basic image of 
the orbit shape that can be recognized in other 
contexts. The drawbacks to this method are that 
it is descriptive rather than quantitative; complex 
orbits usually cannot be simply defined in words; 
and sometimes the adopted terms are not uni- 
versally accepted (e.g., Richstone 1982). Another 
possible problem that arises when dealing with 3D 
orbits is that, while one planar projection may be 
simply described, the other two complementary 
projections may be too complex for words (sec, 
for example. Fig. 9 in Pfenniger 1984). Although 
orbital morphology will not be used as a primary 
characterization tool in this paper, orbit projec- 



tions will be utilized as visual aids when discussing 
most orbits. 

3.3. Spectral Dynamics 

Another way of classifying 2D and 3D orbits is 
the spectral dynamics method first introduced to 
the astrophysical community by Binney & Spergel 
(1982). Briefiy, the spectral dynamics character- 
ization method is based on the fact that regular 
orbits in ncar-intcgrablc potentials exist in re- 
gions of phase space with toroidal topology. As 
noted earlier in footnote 1, here the term "regular" 
means that the number of isolating integrals re- 
spected by the orbit is greater than or equal to the 
dimensionality of the orbital configuration space. 
These phase space tori can be described in terms 
of actions and action angles, so the Fourier spectra 
of the {x{t), y{t), z{t)) components of a regular or- 
bit must consist of discrete lines. The frequencies 
of these lines are linear combinations of the fun- 
damental frequencies associated with the action 
angles. The exact linear combinations of frequen- 
cies contain information about the resonances of 
the orbit. Since ergodic and quasi-ergodic orbits 
do not lie on tori in phase space, their spectra will 
consist of a forest of lines at all frequencies. For 
more complete discussions of spectral dynamics, 
the reader is directed to Binney & Spergel (1982), 
Laskar (1993), Carpintero & Aguilar (1998), and 
Copin et al. (2000). 

Despite the demonstrated utility of the spec- 
tral dynamics method in a variety of contexts, it 
has not been incorporated into the present study 
for two reasons. Primarily, the spectral dynamics 
method does not provide a quantitative measure of 
quasi-ergodic orbits whereas, as mentioned in the 
introduction, realistic galactic potentials probably 
support quasi-ergodic orbits. Secondly, the spec- 
tral dynamics method is strictly applicable only to 
"near-integrable" systems. This means that the 
Hamiltonian for a given system must be express- 
ible in terms of an integrable part plus some small 
perturbation. Numerically determined potentials, 
or their analytical approximations, generally can- 
not be expected to fall into this category. 
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3.4. Lyapunov Exponents 

A more qiiantitative description of an or- 
bit can be made using Lyapunov exponents. 
For a given phase space trajectory, ^{t) = 
(x{t),y(t),z{t),px{t),py{t),pz{t)), the Lyapunov 
exponents measure how a nearby trajectory, 
^{t) + A^(t), where is initially small, diverges 
from the given phase space orbit. A number of pre- 
vious studies of orbits in an astrophysical context 
have utilized Lyapunov exponents to characterize 
the orbits. A few examples that contain excellent 
discussions of the technique are Merritt & Valluri 
(1996), Habib et al. (1997), and especially Udry & 
Pfenniger (1988). Briefly, there is one Lyapunov 
exponent, /cj, for each dimension i of phase space. 
In conservative (Hamiltonian) systems (such as all 
the cases discussed in this paper), these exponents 
are not independent. In fact, J2i = 0- Also, the 
restriction that phase space volumes must remain 
constant (Liouvillc's theorem) connects conjugate 
pairs of exponents. In addition, the exponent that 
corresponds to motion along the direction of mo- 
tion is zero. For example, in a 4D phase space two 
exponents are equal to zero (one is along the di- 
rection of motion and the other is in the conjugate 
direction); the other two exponents are equal in 
magnitude but have opposite signs (Lichtenberg 
& Lieberman 1983). In particular, what is mea- 
sured and referred to as the Lyapunov exponent 
in this paper is the value of the largest positive ex- 
ponent. The attractive attributes of this method 
are its quantitative nature and its applicability to 
orbits in any potential. The main drawback to 
using this method as a solo characterization tool 
is that no distinction can be made between closed 
and unclosed regular orbits. That is, all regu- 
lar orbits display insensitivity to small changes in 
initial conditions and therefore show similar be- 
haviors in their Lyapunov exponents. 

The technique for measuring the largest Lya- 
punov exponent follows from the prescription 
given in Benettin et al. (1976). Lichtenberg & 
Lieberman (1983) give a good introduction to the 
method, starting from the definition of the Lya- 
punov exponent, which is summarized here. The 
assumption is made that two nearby trajectories 
diverge exponentially with time, i.e., 



d{t) = d{0)e 



kt 



(10) 



where d denotes a phase space distance. The Lya- 
punov exponent can then be defined to be. 



k = lim ( — I In 

t->oo,d(0)^0 \ t 



d{t) 



d{0) 



(11) 



If a chosen orbit is quasi-ergodic, so that two 
nearby trajectories diverge as e'^*, then k = a = 
constant. If, instead, the chosen orbit is regular 
and nearby trajectories diverge only as a power 
law in time, d{t) ~ d{0)t", then 



k = - Int. 

t 



(12) 



In this plot of Infc vs. Int should have a 

slope w —1 for large values of t, independent of 
the precise value of the exponent a. 

The definition of k as given in eq.(ll) is com- 
putationally unsatisfactory. Exponential growth 
can quickly lead to numbers that a computer can- 
not represent. Benettin et al. (1976) suggest, in- 
stead, that an orbit be broken into n finite time 
lengths, T (see their Fig. 1 or Fig. 5.6 in Lichten- 
berg & Lieberman 1983 for a pictoral representa- 
tion of this idea). Then, every r time units, the 
distance between neighboring trajectories should 
be re-normalized to the distance between the two 
at the beginning of the orbit. With this technique, 
the Lyapunov exponent is defined to be (Benettin 
et al. 1976), 



1 " 
k„ = — y^ln 

n.T '—^ 



d(iT) 



(13) 



where, in the limit of n — > cxd, fc^ ^ fc. 

The task that remains is to calculate whether or 
not nearby trajectories exponentially diverge. The 
vector of interest here is A^, that is, the vector 
diff'erence between the orbit that is integrated and 
the nearby trajectory. Following Lichtenberg & 
Lieberman (1983), define w = A^. From a linear 
stability analysis (as in Binney & Tremaine 1987, 
§3.5.3), then, 

^=M-w, (14) 

and A4 is a tensor that has components defined 
by, 

dP 

M^-., (15) 

where F = || j , iJ is the system's Hamil- 

tonian, X is the generalized coordinate vector, and 
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p is the conjugate momentum vector. In a sys- 
tem with N degrees of freedom, there are now 3A'' 
equations that must be solved - N for the particle 
trajectory and 2A'' for the phase space difference 
vector. 

The 2 A'' differential equations given by eq.(14) 
are translated into finite-difference equations that 
are solved alongside the equations of motion. For 
orbits that will be analyzed in this paper, fc„ is 
determined for nr = 10^, 5 x 10^, 10^, 5 x 10^, and 
10^. In addition, the value of r has been chosen 
to be the same as At, so the difference vector is 
renormalizcd every timcstcp. The presentation of 
Lyapunov exponents is in the form of plots show- 
ing lnfc„ versus In nr. This form has been chosen 
for these plots because regular orbits (those insen- 
sitive to small changes in initial conditions) have 
slopes ~ — 1, while quasi-ergodic orbits have slopes 
« because they are sensitive to initial conditions. 

3.5. Correlation Integrals 

The main orbital characterization method that 
will be utilized in this paper is the correlation inte- 
gral method. As described by Grassberger & Pro- 
caccia (1983a) and Carnevali & Santangelo (1984), 
the correlation integral provides a measurement of 
phase space dimensionality. The basic idea under- 
lying this method is that for a given orbit, the cor- 
relation integral is calculated for all phase space. 
Then, since a plot of the correlation integral ver- 
sus phase space distance behaves like a power law 
for small phase space distances (Grassberger & 
Procaccia 1983a), the exponent of the power law 
provides a measure of the dimensionality of the 
phase space that is occupied by the orbit. Once 
this dimensionality is known, one can also read- 
ily deduce the number of isolating integrals that 
are respected by the selected orbit. Approaching a 
previously unexamined potential in this way gives, 
for example, basic information about whether or 
not most orbits are regular. This characterization 
then may or may not lead one to perform a more 
complex analysis, such as spectral dynamics. 

Operationally, the steps in determining the di- 
mensionality of a phase space orbit are: 

1. Integrate an orbit for a sufficient number of 
timesteps ("sufficient" will be clarified be- 
low). 



2. Choose a set of sampling points from the or- 
bit. 

3. Calculate the correlation integral, C(r), as 
defined by cq. (16). 

4. Measure a slope from a plot of In C(r) vs. 
Inr. 

Because topics 1 and 2 are somewhat connected, 
they are discussed together. From various trials, 
investigating phase spaces of different sizes, it has 
been determined that approximately 50 orbital pe- 
riods must be completed in order for an accurate 
correlation integral to be found. This means that 
more timesteps are necessary for each 3D orbit 
than for each 2D orbit. With the timestep that 
has been used throughout this investigation, the 
ratio between the number of required timesteps 
is K, 100. (It is likely that these numbers would 
change if a different integration scheme was used.) 
The sampling points are chosen following the sug- 
gestion in Carnevali & Santangelo (1984); specifi- 
cally, random points are selected from a subset of 
the orbital phase space points. The subset con- 
sists of 10^ points taken at equal timestep inter- 
vals. Normally, 10% of these points are randomly 
chosen to be sampling points. While this number 
of sampling points is usually adequate for the pur- 
poses of this study, it generally is best to use as 
many sampling points as is computationally prac- 
tical. 

The correlation integral is determined numeri- 
cally using the following formula (Grassberger & 
Procaccia 1983a): 

N N 

^('■^^iv^^ooiv^^ ^ Q(r--ii:-^-i), (16) 

i=l j=l,^i 

where O is the Heavisidc step function, the are 
phase space position vectors, r is a phase space 
distance, and N is the number of sampling points 
from the orbit. It is now clear why it is advisable 
to use as large a value of N as practical, as the cor- 
relation integral is defined more accurately when 
N is larger. A naive evaluation of eq.(16) is very 
slow for large A^, even if one takes advantage of 
the symmetry of the Heaviside function. Dividing 
phase space into "bins" can significantly speed up 
the evaluation of the double summation. 

The final step in determining the dimensional- 
ity of a phase space orbit is measuring a slope from 
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the plot of In C(r) vs. Inr. A log-log plot is used 
because of the behavior of the correlation integral. 
Grassberger & Procaccia (1983a) state that, for 
r <C 1, C{r) oc r'^ , where v is the dimensionality of 
the phase space orbit. So, when In C(r) is plotted 
against Inr, the slope of any linear section can be 
interpreted as the dimensionality v of the phase 
space orbit within that range of r. In order to as- 
certain the reliability of the values of v that have 
been measured for individual orbits, a minimized 
linear fit to the In C(r) — ln(r) data for five 
independent sets of sampling points has also been 
calculated. Also, at least one order of magnitude 
in r must be covered by the linear section to be 
considered. It is the average slopes {D = (v)) and 
standard deviations (cr) of these five linear fits 
that are reported in the legends of figures such as 
Fig. lib. 

4. Results 

The individual phase space orbits that are dis- 
cussed in this section are listed in Table 1. As 
indicated in the first column of Table 1, the or- 
bits are identified by the potential (or mapping) 
in which they exist. Listed in column 2 are the 
corresponding figure numbers that contain visual 
summaries of the orbital analysis. The initial 
conditions (positions, velocities, and energies) for 
each of the orbits are listed in the next seven 
columns. The last two columns of Table 1 are 
provided as a quick reference to the Lyapunov 
exponent and correlation integral characteriza- 
tions that have been determined for each orbit. If 
the next to last column (labeled Lyapunov) con- 
tains an T (denoting an orbit insensitive to small 
changes in initial conditions), that orbit is regu- 
lar; an 'S' (for sensitive) denotes an orbit that is 
not regular. The last column holds the measured 
dimensionality D of each phase space orbit. 

Before undertaking a discussion of the charac- 
terization of these selected orbits, let us examine 
the expectations for the various models. The 2D 
phase space of the Henon map is not derived from 
a Hamiltonian system, so there are no integrals of 
motion for this system per se. However, each of 
the 4D phase spaces (2D Richstonc, Henon-Heilcs, 
and 2D Gazes bar) is a Hamiltonian system. In 
these cases, the potentials are time-invariant so 
the energy (or ej) must be an integral, and we 



should find that _D < 3 for all orbits. Fully er- 
godic orbits will have D = 3, while quasi-ergodic 
orbits will have 2 < D < 3. Quasi-ergodic orbits 
respect only one full integral of motion, but they 
also have some (unknown) restriction in addition 
to that of their ergodic cousins. Both ergodic and 
quasi-ergodic orbits are sensitive to initial condi- 
tions, so the Lyapunov exponent should remain 
nearly constant with time, as discussed in §3.4. 
On the other hand, a regular orbit must have a 
total number of isolating integrals that is equal to 
or greater than the number of spatial degrees of 
freedom (in this case > 2; specific energy plus 1 or 
2 unknown integrals). So, the allowed phase space 
for each regular (but not periodic) orbit should 
be 2D (4 phase space dimensions minus the 2 iso- 
lating integral dimensions). The behavior of the 
Lyapunov exponent should be that of a regular 
orbit, that is, the slope of lnA;„ vs. Innr should 
be « -1. 

Orbits in 3D Hamiltonian systems occupy 6D 

phase spaces. However, since the potentials are 
still time-invariant, the measured dimensionality 
of the phase space orbit should be, at most, five. 
Regular orbits in 3D potentials must respect at 
least three integrals of motion, D < 3, and as be- 
fore, the Lyapunov exponent should decrease with 
time for a regular orbit. Returning to the nonreg- 
ular orbits, fully ergodic orbits should display a 
dimensionality D = 5, while quasi-ergodic orbits 
should have a noninteger value 4 < _D < 5. Orbits 
with 3 < I? < 4 are neither regular nor quasi- 
ergodic. They respect (at least) two full integrals 
of motion and will be referred to as semiregular 
orbits in this paper. All nonregular orbits have 
Lyapunov exponent plots that have In fc„ vs. In nr 
slopes that are w 0. 

4.1. 2D Phase Space 

As mentioned earlier, the Henon map phase 
space is fractal (see §2.1). The accepted value 
for the dimensionality of this phase space orbit is 
D = 1.25 ±0.02 (Grassberger & Procaccia 1983a). 
Figure 2a shows the {x, y) phase space structure 
of the Henon map; Fig. 2b is a magnified view 
illustrating the fractal nature of this phase space 
orbit. Figure 2c shows the In C(r) vs. Inr plot 
(hereafter, the C — r plot) that was obtained in 
this study for the Henon map using the correlation 
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integral method. The resuh is D = 1.21 ± 0.01. 
This is also the value that Grassberger & Procac- 
cia (1983a) obtained using the correlation integral 
method. The difference between this result and 
the accepted value is discussed in Grassberger & 
Procaccia (1983b). They claim that there is a sys- 
tematic error in the determination of the Hcnon 
map correlation integral. Interestingly, this error 
does not occur in correlation integral determina- 
tions for other maps such as the Kaplan- Yorke 
map, the logistic equation, or the Lorenz equation 
(Grassberger & Procaccia 1983b). They conclude 
that the error in the Hcnon map result arises 
from the sample phase space points not uniformly 
covering the phase space orbit. However, their 
suggested remedy does not work well when phase 
space is populated by numerically integrated or- 
bits. In the following sections it will be shown 
that this slight disagreement does not diminish 
the quantitative usefulness of the correlation inte- 
gral method as a tool for characterizing orbits. 

4.2. 4D Phase Space 

For orbits in 2D potentials, the presentation 

of results will adhere to the following form: In 
each figure, the frame labeled (a) contains the or- 
bital trajectory; frame (b) displays the C — r plot 
obtained using the correlation integral method; 
frame (c) shows the Lyapunov exponent plot as 
lnfc„ vs. In nr; and frame (d) contains either the 
{x — px) or [R — pji) surface of section diagram 
for the orbit. The legends of frames labeled (b) 
also will provide quantitative information derived 
from the C — r plot as discussed in §3.5, namely 
the slope D and standard deviation a. In most 
Lyapunov exponent plots, a dot-dashed line with 
a slope of -1 also is included because, as discussed 
in §3.4, the behavior of kn should approach this 
slope as nr —> 00 if an orbit is regular. 

The results for four orbits in the 2D Richstone 
potential are shown in Figs. 3 through 6. In ev- 
ery case, all three techniques for characterizing 
the orbits indicate that the orbits are regular: the 
surface of section plots are composed of invariant 
curves; the Lyapunov exponent drops with time 
with a slope of minus one; and from the corre- 
lation integral, the measured dimensionality has 
a value < 2. Note that the Lyapunov exponent 
plots in Figs. 3c and 6c are basically identical and 



therefore make no distinction between the closed 
and unclosed regular orbits. However, the differ- 
ence is clearly illustrated by the differing slopes of 
the C — r plots in Figs. 3b and 6b: the periodic 
orbit displays D = 1, instead of D = 2. 

At this point, it is worth noting some general 
features of C — r plots derived from orbits (of any 
dimensionality) as opposed to mappings. As the 
value of r increases, the plotted points will often 
deviate from the line of slope D. This is because 
a clean linear relation is expected only for r <C 1. 
Another way to think about this is that, for a 
sufficiently large value of phase space distance 
ro, the entire phase space orbit will be enclosed. 
Then, for r > ro, C{r) = 1.0, so the C - r plot 
must tend toward a slope of zero at large enough 
values of r. Deviations from a linear slope of D 
may also occur at the smallest values of r, but 
for a different reason. Since only a finite number 
of sampling points is used, there is necessarily a 
lower limit to the smallest distance that can be 
measured between any two points. Basically, this 
is a small-number statistics problem. As a larger 
number of sampling points is used, the linear fit 
generally becomes tighter and extends to smaller 
values of r. 

Four orbits with e = 1/8 in the Henon-Heiles 

potential are displayed and analyzed in Figs. 7 
through 10. Unlike the 2D Richstone potential, 
the Henon-Heiles potential supports some quasi- 
ergodic orbits at this energy. All three methods 
of characterization indicate that the orbits shown 
in Figs. 7 and 8 are regular, while the ones shown 
in Figs. 9 and 10 are quasi-ergodic. Focusing on 
the two quasi-ergodic orbits, notice that, (1) the 
Lyapunov exponent is approximately constant, 
signaling an exponential departure of two initially 
neighboring trajectories; (2) the surface of sec- 
tion is no longer composed of a smooth invariant 
curve; and (3) the C — r plot identifies a dimen- 
sion greater than two. We know that the orbits 
shown in Figs. 9 and 10 are quasi-ergodic rather 
than fully ergodic because of the measured dimen- 
sionality. This is a distinction that can be made 
clearly from the C — r plot, but it is not possible 
from a measurement of the Lyapunov exponent 
alone. The deviations from D at small values of r 
in Figs. 9b and 10b are simply artifacts of the lim- 
ited number of orbital timesteps. This is similar 
to the small number statistics problem that was 
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discussed above. When the number of timesteps 
is increased, these deviations disappear and the 
line with slope D extends to smaller r values. 

Four different orbits that arc supported by the 
2D analytical Gazes bar potential are presented 
here in Figs. 11 through 14. Once again, all 
three methods of characterizing these 2D orbits 
agree: the orbits shown in Figs. 11, 12, and 13 are 
regular, while the one shown in Fig. 14 is quasi- 
ergodic. The orbit shown in Fig. 11a has £) ?a 1, 
rather than D = 2; hence, it is periodic. This orbit 
is the parent of what has been called the 'bow tie' 
family of orbits in Barnes & Tohline (2001). The 
orbit shown in Fig. 12a appears to be trapped near 
the periodic bow tie orbit. The most interesting 
aspect of this orbit is visible in the C — r plot, Fig. 
12b. There are two linear sections in the C — r 
plot. The one that exists for small r has a slope of 
D K 2. At larger r, however, the linear section has 
a slope oi D Ki 1. This slope discontinuity persists 
even when this orbit is followed through 100 times 
as many timesteps. This suggests that the dis- 
continuity and the r location of the discontinuity 
are physically relevant (as opposed to the cases in 
the previous paragraph where differing slopes were 
numerical artifacts). We interpret this difference 
using the following analogy presented by Gleick 
(1987). Imagine viewing a ball of string from very 
far away. If asked to describe the dimensionality 
of the ball of string, which appears to be a point, 
the answer would be zero. Moving closer, the ball 
can be seen to be extended. The ball now appears 
to have dimensionality equal to 2. Moving even 
closer, the ball now appears to have a finite extent 
in a direction perpendicular to the dimensions al- 
ready present as well, and is thus a 3D object. 
However, observing the ball of string at very close 
range, the one-dimensional nature of the string 
becomes apparent. So, the dimensionality of the 
phase space orbit can change depending on what 
length scale is observed. Indeed, for small enough 
values of r, all phase space orbits derived from 
particle trajectories are really ID objects. With 
this in mind, the segment with slope D 1 \s 
interpreted as demonstrating that the orbit is not 
far from being closed. The r value at which the 
change in slope occurs identifies a critical length 
scale for this orbit that cannot be determined via 
the Lyapunov exponent method or from surfaces 
of section. 



The quasi-ergodic orbit shown in Fig. 14 has a 
dimensionality D k. 2.54. Unlike the discontinuity 
present in the regular orbit C — r plot (Fig. 12b), 
the apparently linear slope at smaller r values is 
dependent on the number of timesteps taken for 
the orbit, as with the orbits shown in Figs. 9 and 
10. 

4.3. 6D Phase Space 

As with the 2D orbits previously discussed, the 
figures containing results for individual 3D orbits 
all have a similar form. Each figure consists of: 
a frame labeled (a) illustrating the x — y projec- 
tion of the orbit; a frame labeled (b) illustrating 
the X — z orbital projection; a frame labeled (c) 
illustrating the y — z projection of the orbit; a 
frame (d) showing the C — r plot for the orbit; 
and finally, a frame labeled (e) that displays the 
Lyapunov exponent plot derived from the orbit. 

Figures 15 through 17 display results for three 
different orbits in the 3D Richstone potential. As 
expected, these results support the characteri- 
zation of these orbits as regular. In this case, 
the three integrals of motion are the energy, the 
z-componcnt of angular momentum (which was 
explicitly conserved in the 2D case), and the un- 
known integral that was also apparent in the 2D 
case. These figures demonstrate that the correla- 
tion integral method agrees with other character- 
ization techniques for 3D as well as 2D orbits. 

While it has been important to use individ- 
ual orbits to compare characterization methods, 
it is useful to investigate ensembles of orbits us- 
ing only correlation integrals. As stated in §2.4, 
Stackel potentials provide an excellent test bed 
for characterization methods. One hundred orbits 
of varying energies (25 orbits each at E = —0.2, 
—0.4, —0.6, —0.8) have been integrated in the po- 
tential given by eq. (7). These orbits have zero 
initial velocities and initial positions that cover 
one octant of the given energy surface. 

The results are shown in Fig. 18. This his- 
togram shows the number of orbits versus the 
number of integrals of motion respected by those 
orbits. Note that the number of integrals of mo- 
tion / is determined by subtracting the D value 
taken from a C — r plot from the phase space 
dimension. In this case, I = Q — D. There is 
some scatter about the expected value of / = 3. 
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This scatter is due to the fact that only 10,000 
samphng points have been used. For example, the 
orbit that here appears to respect only 2.6 inte- 
grals of motion gives a value of 3 when 10 times 
as many sampling points are used. However, this 
improved accuracy comes at the price of comput- 
ing time. Ten times as many sampling points 
means an increase by a factor of 100 in computing 
time (from about 7 minutes to about 12 hours per 
orbit) . This made it impractical to use more sam- 
pling points for each of the 100 orbits. Despite 
this scatter, the majority (81%) of the results lie 
between 2.8 and 3.1. The small cluster of orbits 
that respect 4 integrals of motion are those con- 
fined to symmetry planes of the potential. 

Another histogram of an ensemble of 100 orbits 
is shown in Fig. 19. These orbits are supported by 
the 3D analytical Gazes bar potential and, like or- 
bits in the Stackel potential, there are four groups 
of 25 orbits at different energies (ej = —0.96, 
-0.85, -0.75, -0.63). The initial positions of 
these orbits have been chosen to cover one quarter 
of the given energy surface. The initial velocities 
of these orbits have been specified by the Restric- 
tion Hypothesis discussed in §3.1. 

The most conspicuous feature of this histogram 
is that the integrals of motion are not concentrated 
around a single value. What this histogram shows 
is that the 3D Gazes bar potential supports or- 
bits that respect 3 integrals of motion (regular), 
2 integrals of motion (scmircgular) . and 1 integral 
of motion (quasi- ergo die). While this histogram 
has not been created with a self-consistent stellar 
distribution (which is a future project), it does 
give the flavor of how such a distribution will ap- 
pear. Additionally, histograms of the 25 orbits at 
the various Jacobi constants are also interesting. 
They show a trend for an increasing fraction of 
nonregular orbits with increasing energy. This is 
similar to the trend for orbits in the Henon-Heiles 
potential. 

5. Discussion & Conclusions 

Building on the work of Grassberger & Procac- 
cia (1983a) and Garnevali & Santangelo (1984), 
a flexible technique for characterizing orbits in 
3D potentials known as the correlation integral 
method has been presented. This method ana- 
lyzes phase space orbits and returns a single num- 



ber, the dimensionality of the phase space orbit. 
From this number and the dimensionality of the 
underlying phase space, the number of isolating 
integrals of motion respected by an orbit can be 
determined. This number can then be used as a 
quantitative characterization attribute. 

The implementation of the correlation integral 
method has been tested for a variety of previously 
studied systems. More familiar characterization 
tools, such as surfaces of section and Lyapunov 
exponents, support the results obtained with the 
correlation integral method. The advantages of 
the correlation integral are most apparent when 
used to characterize orbits in 3D potentials. A 
unique, 3D potential based on a numerically cre- 
ated potential-density pair called the Gazes bar 
has been presented. The Gazes bar is rotating and 
has no geometrical symmetries that give rise to an- 
alytical integrals of motion. However, the corre- 
lation integral method demonstrates that regular 
orbits exist in the Gazes bar potential. Addition- 
ally, the correlation integral method distinguishes 
between orbits that respect two integrals of mo- 
tion and those that respect only one integral. 

The simple fact that the correlation integral 
method can reproduce the results of other char- 
acterization methods is not enough to warrant its 
adoption. Here the various categorizing methods 
discussed in this paper are compared and con- 
trasted. 

• When analyzing orbits in a 4D phase space, 
surface of section diagrams are the sim- 
plest and clearest way to characterize or- 
bits. However, there is no simple quanti- 
tative measure that describes quasi-ergodic 
orbits in surface of section diagrams. Also, 
these diagrams arc not easily translated to 
6D phase spaces. The correlation integral 
method addresses both of these problems. 

• Lyapunov exponents provide quantitative 
measures of orbital regularity in arbitrary 
2D and 3D potentials. Unfortunately, all 
regular orbits, closed and unclosed, share 
the same signature in Lyapunov exponent 
plots. For orbits in 3D potentials, the be- 
havior of the Lyapunov exponent is also the 
same for all nonregular orbits. That is, no 
distinction is made between orbits that re- 
spect only one integral of motion and those 
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that respect two. The correlation integral 
method distinguishes between these types of 
regular (periodic and unclosed) and nonreg- 
ular (respecting only one or two integrals of 
motion) orbits. 

• The spectral dynamics characterization 
method rests on the assumption that the 
potential being investigated is at least near- 
integrable. When this is true, the spectral 
dynamics method provides an excellent way 
to analyze orbits. If, instead, the poten- 
tial is the result of a numerical simulation 
or is not clearly near-integrablc, then spec- 
tral dynamics may not give relevant results. 
Also, spectral dynamics only returns infor- 
mation about regular orbits. If it is true that 
quasi-ergodic orbits play a role in galactic 
structure, then it will be necessary to utilize 
a characterization tool, such as the corre- 
lation integral method, that incorporates 
those orbits. 

The correlation integral method could prove 
useful for a variety of future studies, but there 
are two particular avenues of interest. Investiga- 
tions of the impact of massive central objects on 
stellar orbits could benefit from the correlation 
integral method in the following way. By creating 
histograms, like those shown here in Figs. 18 and 
19, from orbits in potentials with point masses, 
the general characteristics of such a stellar system 
can be seen. This information may prove useful 
when creating self-consistent stellar distribution 
functions to describe black hole/steUar systems 
(e.g., van der Marel et al. 1998, Gebhardt et al. 
2000b). 

The correlation integral method also should be 
a good way to analyze stellar orbits in cosmolog- 
ical simulations that include star formation. By 
simply keeping a record of the phase space coor- 
dinates of all stars during a simulation, the corre- 
lation integral method can be used to analyze the 
stellar distribution functions in the emerging po- 
tentials. To use the Lyapunov exponent method, 
extra equations would need to be solved during an 
already complex calculation. Similarly, the spec- 
tral dynamics technique may not be expedient if 
the potentials that arise are not sufficiently close 
to integrable forms. 

In conclusion, the correlation integral method 



is a flexible, quantitative, and straightforward way 
of characterizing orbits in 3D potentials. It is rec- 
ommended that it be broadly adopted as a tool 
for characterizing the properties of orbits in all 
Hamiltonian dynamical systems. However, there 
are three specific cases of astrophysical interest 
to which the correlation integral method seems 
particularly well suited: analyzing stellar distri- 
bution functions in analytically and numerically 
spccifled models of steady-state galactic potentials 
(especially those with Hamiltonian chaos); investi- 
gating the orbital structure supported by galactic 
potentials formed in cosmological simulations; and 
quantifying the response of stellar systems to po- 
tentials that contain central point masses. 
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"^The letters that appear in this column distinguish between orbits that are insensitive (I) to small changes 
in initial conditions and those that are sensitive (S) to such changes. All orbits labeled T are regular. 
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Fig. 2. — Henon map. (a) Plot of phase space 
points obtained after 100,000 iterations of eqs. (2) 
and (3). (b) Magnification of the plot in (a). The 
small scale substructure is a telltale sign of the 
fractal nature of the phase space orbit, (c) A plot 
of the correlation integral (In C(r)) versus phase 
space distance (Inr) for this mapping. The value 
of D shown in the legend has been measured from 
the slope of the line; a is the standard deviation 
measured for D, as described in §3.5. 



Fig. 3. — Orbit #1 in the 2D Richstone potential; 
see Table 1 for initial conditions, (a) The R — z 
trajectory of this orbit, (b) The associated C — r 
plot listing the average slope (D) and standard 
deviation (cr). (c) The plot of the Lyapunov ex- 
ponent (lnfc„) versus integration time (In nr) for 
this orbit, (d) The {R — pn) surface of section for 
this orbit. 



Fig. 4. — Orbit #2 in the 2D Richstone potential; 
each frame contains information as described in 
the caption to Fig. 3. 



This 2-column preprint was prepared with the AAS lAlJjjX 
macros v5.0. 



Fig. 1. — (a) Equipotential contours of the nu- 
merical Gazes bar in the x — z plane for a; > 0; 
(b) equipotential contours of the analytical Gazes 
bar in the same plane, (c) Equipotential contours 
of the numerical Gazes bar in the y — z plane for 
2/ > 0; (d) equipotential contours of the analytical 
Gazes bar in the same plane, (e) Equipotential 
contours of the numerical Gazes bar in the x — y 
plane; (f) equipotential contours of the analytical 
Gazes bar in the same plane. 



Fig. 5. — Orbit #3 in the 2D Richstone potential; 
each frame contains information as described in 
the caption to Fig. 3. 



Fig. 6. — Orbit #4 in the 2D Richstone potential; 
each frame contains information as described in 
the caption to Fig. 3. 
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Fig. 7. — Orbit #1 in the Henon-Heiles potential; 
see Table 1 for initial conditions, (a) The R — z 

trajectory of this orbit, (b) The associated C — r 
plot listing the average slope {D) and standard 
deviation (a), (c) The plot of the Lyapunov ex- 
ponent (lnfc„) versus integration time (Innr) for 
this orbit, (d) The {R — pn) surface of section for 
this orbit. 



Fig. 8. — Orbit #2 in the Henon-Heiles potential; 
each frame contains information as described in 
the caption to Fig. 7. 



Fig. 15. — Orbit #1 in the 3D Richstonc potential; 
see Table 1 for initial conditions, (a) The x — y 
plane projection of this orbit, (b) The x — z plane 
projection of the same orbit, (c) The y — z plane 
projection of the same orbit, (d) The associated 
C — r plot listing the average slope {D) and stan- 
dard deviation {a), (e) The plot of the Lyapunov 
exponent (InA;^) versus integration time (Innr) for 
this orbit. 



Fig. 9. — Orbit #3 in the Henon-Heiles potential; 
each frame contains information as described in 
the caption to Fig. 7. 



Fig. 10. — Orbit #4 in the Henon-Heiles potential; 
each frame contains information as described in 
the caption to Fig. 7. 



Fig. 16. — Orbit ^2 in the 3D Richstonc potential; 
each frame contains information as described in 
the caption to Fig. 15. 



Fig. 11. — Orbit #1 in the 2D Gazes bar potential; 
see Table 1 for initial conditions, (a) The x — y 

trajectory of this orbit, (b) The associated C — r 
plot listing the average slope {D) and standard 
deviation (ct). (c) The plot of the Lyapunov ex- 
ponent (lnfc„) versus integration time (Innr) for 
this orbit, (d) The {x — Px) surface of section for 
this orbit. 



Fig. 17. — Orbit in the 3D Richstonc potential; 
each frame contains information as described in 
the caption to Fig. 15. 



Fig. 12. — Orbit #2 in the 2D Gazes bar potential; 
each frame contains information as described in 
the caption to Fig. 11. 



Fig. 18. Histogram of the number of orbits ver- 
sus the number of integrals of motion for 100 orbits 
in the Stackel potential described in §2.4. 



Fig. 13. — Orbit ^^3 in the 2D Gazes bar potential; 
each frame contains information as described in 
the caption to Fig. 11. 



Fig. 14. — Orbit #4 in the 2D Gazes bar potential; 
each frame contains information as described in 
the caption to Fig. 11. 



Fig. 19. Histogram of the number of orbits ver- 
sus the number of integrals of motion for 100 orbits 
in the 3D analytical Gazes bar potential. 
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